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The momentum difTusion coefficient for heavy quarks is studied in a deconfined gluon plasma 
in the static approximation by investigating a correlation function of the color electric field using 
Monte Carlo techniques. The diffusion coefficient is extracted from the long distance behavior of 
' such a correlator. For temperatures Tc < T ^ 2Tc, our nonperturbative estimate of the diffusion 

^Sj , coefficient is found to be very different from the leading order perturbation theory, and is in the 

^ ■ right ballpark to explain the heavy quark fiow seen by PHENIX at RHIC. 
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I. INTRODUCTION 



The charm and the bottom quarks arc very important tools in our quest to understand the nature of the quark- 
Oh' gluon plasma created in the relativistic heavy ion collision experiments. Since the masses of both of them are much 
^ ■ larger than the temperatures attained in RHIC, and in LHC, one expects these quarks to be produced largely in the 
early pre-equilibrated state of the collision, and thus provide a window to look into the early stages of the fireball. 
Furthermore, perturbative arguments suggest that the energy loss mechanism for energetic heavy quarks in medium 
should be different from that of the light quarks. A comparative study of the energy loss for the heavy and light 
. quark jets therefore leads to crucial insights into the way the quark-giuon plasma interacts. 

' For light quark jets, gluon radiation ( "bremsstrahlung" ) is expected to be the leading mechanism for energy loss in 
J — [ medium It has been argued that gluon bremsstrahlung is suppressed for jets of heavy quarks 0], and collisional 
■ energy loss may be the dominant mechanism for thermalization of not-too-energetic heavy quark jets [1, Since 



collision with a thermal quark does not change the energy of a heavy quark substantially, one would expect that the 
. thermalization time of the heavy quarks is much larger than that of the light quarks. As most of the elliptic flow is 
developed early, the azimuthal anisotropy parameter, V2, of the hadrons with heavy quarks can be expected to be 
much less than that of the light hadrons. 

Interesting predictions follow from these simple, weak coupling-based intuitions, which can be, and have been, 
checked in the RHIC experiments. One expects a mass ordering of the elliptic flow: V2 ^ V2 ^ V2 . Here 
h, D, B refer to the light hadrons, mesons of the D family (one charm and one light quark) and those in the B 
\ family (one bottom and one light quark). The nuclear suppression factor, RaAi is also expected to show a hierarchy: 
^'aa ^ ^ ^AA- Experimentally, on the other hand, it was found that the heavy flavor mesons show a large 

elliptic flow, V2 ^ i'2i ^'^d a strong nuclear suppression, i?^^ ^ ^\a^ nuclear suppression being comparable to 
that of tt" for pT > 2 GeV 0,0. 

Even if the kinetic energy of the heavy quark is 0(T), where T is the temperature of the fireball, its momentum will 
be much larger than the temperature. It is, therefore, changed very little in a single collision, and successive collisions 
can be treated as uncorrelated. Based on this picture, a Langevin description of the motion of the heavy quark in the 
medium has been proposed 0, 0, 01- ''^z, the elliptic flow parameter, can then be calculated in terms of the diffusion 
coefflcient of the heavy quark in the medium. The diffusion coefficient has been calculated in perturbation theory 
[li 0| • While the experimental results for the elliptic flow of the charmed mesons and its pr dependence seem to be 
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well described by this formalism for moderate pt ^ 2 GeV, the value of the diffusion coefficient needed to explain 
the experimental data is found to be at least an order of magnitude lower than the leading order (LO) perturbation 
theory (PT) result [131 Furthermore, the LOPT value itself leads to heavy flavor flow which falls way short of the 
data. The contribution of the next-to-leading order (NLO) in perturbation theory has been calculated recently [^j- 
Although, it was found to change the LO result by a large factor at temperatures < 2Tc, this should perhaps be taken 
as an indication of the inadequacy of perturbation theory in obtaining a reliable estimate for the diffusion coefficient 
in the temperature range of interest. 

A nonpcrturbative estimate of the diffusion coefficient, D, in QCD is, therefore, essential to understand the heavy 
quark flow in the Langevin formalism. Lattice QCD, together with numerical Monte Carlo techniques, provides 
the only way of doing first principle nonperturbative calculations in the quark-gluon plasma. Unfortunately, such 
calculations are done in Euclidean space, and extracting a real time object like the diffusion coefficient requires an 
analytic continuation, which is extremely difficult. However, estimation of various transport coefficients have already 
been attempted, with varying degrees of success [lo| . In order to estimate the heavy quark diffusion coefficient, 
one can study the correlator of the heavy quark current, Qj^Q. Early attempts to extract D this way showed that 
the correlator has very little sensitivity to D Preliminary results of a calculation of D extracted from Qj'^Q 

correlator, using much finer lattices than used before, have been presented recently in the temperature range between 
1.5 — 3 Tc for a gluon plasma [l^. They do find a value which is much lower than the LOPT result, and in the right 
ballpark to explain the experimental heavy quark flow. 

Two of the difficulties in extracting the diffusion coefficient from Qj^Q correlator are: i) the behavior of the structure 
of the spectral function near the lj ^ 2toq regime can affect the structure at low w, and ii) the diffusion coefficient is 
obtained from the width of the narrow transport peak at w — > 0, which is difficult to extract. In the infinitely heavy 
quark limit, another approach to the diffusion coefficient has been suggested in Refs. 0, In this static limit, the 
propagation of heavy quarks is replaced by Wilson lines, and the formalism of Ref. [3| reduces to the evaluation of 
retarded correlator of electric fields connected by Wilson lines [l^ . In Ref. [I3|] , this formalism was used to calculate 
the diffusion coefficient for the AA = 4, SU(A^c — > oo) gauge theory, using the AdS/CFT correspondence. A parametric 
dependence on coupling very different from weak coupling perturbation theory was obtained. On the other hand, for 
the pure SU(3) gauge theory, a leading order perturbative calculation of this correlator led to a negative value for the 
diffusion coefficient at moderate temperatures I la. 



The formalism outlined in Ref. [ij] is suitable for Monte Carlo calculation on the lattice. As we outline in the next 
section, this involves the calculation of Matsubara correlators of color electric field operators, and extracting the low 
frequency part of the spectral function from it. Of course, the usual problems of extraction of the spectral function 
from the Matsubara correlator mean that calculation of the diffusion coefficient remains a highly nontrivial task. A 
first attempt to calculate the diffusion coefficient from the electric field correlator lead to very large values of the 
diffusion coefficient, close to the perturbation theory value (l6j . Preliminary results from a recent calculation [l7j , on 
the other hand, gave values in the temperature range 1.5-3 Tc close to the experimental results. 

In this work, we use the formalism of [ij, [l3l to calculate the diffusion coefficient of the deconfined gluonic plasma 
in a moderate temperature range, Tc < T ^ '^Tc. The aim is to understand whether a small diffusion coefficient, as 
found in the analysis of the experimental data is consistent with QCD. The plan of the paper is as follows. In 
the next section we outline the formalism. In Sec. IIIII we explain the operators and the algorithm. Sec. IIVI has our 
results. A discussion of the results, including their connection with experiments, is contained in Sec. [V] Some details 
of Sees. Hn] and |IV] are relegated to the appendices. 



In this section, we outline the formalism of Refs. [13|, |14[ in more detail. We first sketch the arguments leading 
to the Langevin formalism, and then discuss the quantum field theoretic calculation of suitable correlation functions. 
This discussion closely follows Refs. [1, Then we discuss the issues related to the extraction of the diffusion 

coefficient from the Matsubara correlator. 

It is easy to see why the motion of a quark much heavier than the system temperature can be described in the 
Langevin formalism. If the kinetic energy is ~ T, then the momentum, ^ \J MT , is not changed substantially in 
individual collisions with thermal giuons and quarks, which can only lead to a momentum transfer ^ T . Therefore, 
the motion of the heavy quark is similar to a Brownian motion, and the force on it can be written as the sum of a 
drag term and a "white noise" , corresponding to uncorrelated random collisions: 




II. 



FORMALISM 



dt 



VDPi + 



K dij 5{t~t'). 



(1) 



3 



From Eq. ([T]) the momentum diffusion coefBcient, k, can be obtained from the correlation of the force term: 



The drag coefficient, ry^), can be connected to the diffusion coefficient using standard fluctuation-dissipation relations 



Here M is the heavy quark mass. 

To have a field theoretic generalization of Eq. ^ , one flrst introduces the conserved current for the heavy quark 
number density, J^{x^t) = ip{x,t)j'^'ip{x,t), where '0 is the heavy quark field operator. The force acting on the heavy 
quark is given by M dJ^/dt and so, Eq. ^ generalizes to 



— lim 

3 w^o 



M2 

lim — — — 



j(t-t') 



d^a 



1 [ d.r{x,t) d.r{Q,t') 

2 I Jt ' dt' 



(4) 



where x is the spatial integral of the density correlator: 



00 



d^x{J°{x,t) J°(0,t)) = Tx 



and is directly proportional to the number density for a system of non-relativistic quarks. 

Since we are working in the heavy quark limit, the force term and the number density term are easy to infer: 



(5) 



= {</.t£;v - 0^E'0}, 

J° = (t)U + S^S (6) 

where (p and 9 are the two-component heavy quark and antiquark field operators, respectively, and is the color 
electric field. In leading order expansion in 1/M, only the electric field contributes to the force term. 

With the substitution of Eq. (jH]), the real time correlator in Eq. (U) can be calculated as the analytical continuation 
of the Matsubara correlator, 

Ge{t) = -^E ,1™, ^ / '^'^ - O^E'e){T,x) _ e^E'e){{),Q))- (7) 

The spectral function, p{uj), for the force term is connected to Ge{t) by the integral equation [l^ 



Ge{t) = / — p{uj) ^. , ■ (8) 



TT sinh 



2T 

The momentum diffusion coefficient, Eq. ([4]), is then given by 



2T 

K = lim — p{uj)- (9) 
tj— >o uj 

Since we are working in the limit of infinitely heavy quarks, the expression ([7]) simplifies considerably. The heavy 
quark correlators give a static color field, besides an exponential suppression factor coming from the heavy quark 
mass: {9a{T, x)9l{0,0)) = d^{x) Uab{T,0) exp(— Mr), where Uab{T,0) is the timelike gauge connection, and the delta 
function comes because the infinitely heavy quark does not move spatially. The exponential factor cancels with a 
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similar factor from x i resulting in a rather simple expression for the infinitely heavy quarks: 



G^"*(t) = E (l^c^'l \uW,r) Mt,0) C/(t,0) E,iO,0) 



(10) 



where L = tr U{P, 0) is the Polyakov loop. Once again, intuitively it is easy to understand Eq. pU)) : for the infinitely 
heavy quarks, all that the force-force correlator gives is the correlator of color electric fields, connected through Wilson 
lines, and normalized by the Polyakov loop. 

In order to connect G]f^{T) measured on the lattice to physical correlator of electric fields, wc need to multiply by 
a renormalization factor: 



GE{r) = Z{a)G%-\T) (11) 

where Z{a) = Z'^ is the lattice spacing dependent renormalization factor for the electric field correlator. A nonper- 
turbative evaluation of the electric field operator used here is not available. However, the renormalization factor is 
expected to be dominated by the self energy correction, which can be taken into account by a tadpole correction [l^ . 
In fact, with other discretizations of the electric field operator it has been found that the tadpole factor gives a very 
close approximation to the nonperturbative renormalization factor [20| . Here we use the tadpole factor to renormalize 
the electric field. 

The extraction of p{uj) from Ge{t) using Eq. ([8|) is an extremely difficult problem. In general, the kernel in Eq. 
[5] will have zero modes on a discrete lattice, making the inversion problem ill-defined. In the ideal case, one may be 
able to impose some reasonably general conditions on p(a;) and be able to make the problem invertible. However, in 
situations like ours, where Ge{t) is measured only on 0(10) data points with errors, the problem of extraction of 
p{ui) becomes a completely ill-posed problem without any further input. 

For some problems, a Bayesian analysis, with prior information in the form of perturbative results, has been useful. 
In general, though, a stable application of these techniques require both a very large number of points in the r 
direction and very accurate data for Ge{t). For the kind of extended objects we are considering, wrapping the lattice 
in the Euclidean time direction, it is very difficult to obtain both together, as the error on the correlators grows with 
the number of points in the r direction. 

Parameterizing p{uj) in terms of a small number of parameters, therefore, seems to be a simple way to make the 
inversion problem well-posed. In our case, the leading order perturbative form of the spectral function ^ 6a;^. Also in 
the oj — >■ regime, we need p{u)) nw to get a physical value of the diffusion constant using Eq. The calculation 
of Ref. [l^ got p{uj) = cuj for the M — A supersymmetric Yang-Mills theory. Motivated by this, we use a simple 
ansatz for the spectral function, 

piiuj) = aa;e(w-A) + buj^ (12) 

We do not include any running in b, which is proportional to as in the leading order. This approach is similar in 
spirit to that used in Ref. [2l| to calculate electrical conductivity. Note that the NLO PT calculation of Eq. ([TU]) 
leads to a negative value of k, and in general, seems to deviate more from the lattice correlators than the LO result. 
So we use the LO form for the high uj part. 

On the other hand, in the calculation with classical lattice gauge theory the spectral function of the color 
electric field was found to have the behavior 



p{uj) ^ ctanh-^ for uja<^l. 
So to crosscheck the dependence on our assumption, we also use a second fit form, 

P2{uj) = ctanh^^ e(w - A) + bu^. (13) 

In practice, we use these postulated forms for p{ll!) to evaluate Ge{t), and fit it to the long distance correlation 
function measured on the lattice. At large w, of course, this form is not valid, and a complicated form, that takes into 
account the effect of the lattice Brillouin zones, will have to be considered. We tried using the free lattice spectral 
function instead of the term in Eqn. p2|). However, that did not improve the fit quality and in particular, did not 
seem to capture the very short distance behavior of the data any better. So in this work, we restrict ourselves to the 
large distance regime in our fits, and expect that in this regime our simple form will suffice for a first estimate of the 
diffusion coefficient. 
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FIG. 1: Illustrating the use of the multi-level algorithm for the calculation of G£;,num(4) on a A^t = 12 lattice. 



III. NUMERICAL DETAILS 



For the lattice evaluation of the correlator Ge{t), 
Following Ref. [3l we choose the discretization 



we first need to choose a discretization of the electric field. 



E,{x,t) = U,{x,t) Ui{x + i,T) - Ui{x,T) U,{x + l) 

which is a direct latticization of the relation Ei = [Dq^ Di]. As Ref. [3l suggests, this form of the discretization of the 
electric field is expected to be less ultraviolet sensitive than the more common discretization in terms of the plaquette 
variable. 

The numerator of Eq. (fTO]) can then be written as 



t-l t+T-l P-I 

C'{t) = II Ui{xi) ■ U,{t) ■ J] t/4(x4) • U}{t + T) ■ J] Ui{xi) 



(14) 



X4 — t 



X4—t-\-T 



The evaluation of C*(t), Eq. (|T4|). is known to be difficult for large r, because the signal-to-noise ratio decays 
exponentially. The multilevel algorithm [2^ was indeed devised to take care of such problems. We adapted it for 
calculation of the electric field correlation functions. The lattice is divided into several sublatticcs. The expectation 
value of the correlation functions are first calculated in each sublattice by averaging over a large number of sweeps in 
that sublattice while keeping the boundary fixed. A single measurement is obtained by multiplying the intermediate 
expectation values appropriately. The number of sublatticcs and the number of sublattice averagings were tuned for 
the various sets, so as to get correlators with % level accuracy. An explicit example is shown in Fig[l]which illustrates 
the calculation of num(4) on a iVt = 12 lattice with four sub-lattices, each with a thickness of three lattice spacings. 
It is important to note that one needs to store all the intermediate sub-lattice averages separately before they can 
be multiplied at end of the update of the whole lattice to construct the correlation functions. This imposes memory 
constraints for simulating large lattices. 

The advantage of the multilevel algorithm can be seen from the following estimate: for /? = 6.9, A''^ = 20 and 
Ng = 36, the correlator for t = 3a, G£;(3) has the value of 1.317(2) xlO~^ from 350 multilevel measurements. 
The multilevel algorithm takes about 800 minutes to yield a single measurement on an Intel Xeon CPU processor 
with a speed of 2.5 GHz. For the same correlator, the standard method, using an updating with a combination of 
overrelaxation and heatbath steps, led to a value 1.2(2) xlO~^ for a runtime of about 8500 minutes on the same 
machine. Using the usual 1/Vi dependence of the error on nmtime, the multilevel algorithm is seen to be about 300 
times more efficient than the standard algorithm for G£;(3) for this lattice. The efficiency of the multilevel algorithm 
increases significantly for larger values of t. A similar comparison for G£;(10) gives a factor of about 2000 (order of 
magnitude larger) relative efficiency for the multilevel algorithm. Thus, use of the multilevel scheme is indispensable 
for calculations at the larger values of r 31| , since these are required to be known with high precision for the extraction 
of the diffusion coefficient. 

To get the results for various temperatures and volumes, we ran our simulations at a number of bare couplings with 
A't = 12 - 24 and Ng/Nt = 2-4, for temperatures from just above Tc to 3 Tc- A reliable extraction of the diffusion 
coefficient was possible, however, only for lattices with Nr > 20. A list of such lattices used by us is given in Table 
|T] below. To obtain the temperature scale, we follow the strategy outlined in We calculate at each /3 from 
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FIG. 2: (Left) Ge{t) for one of our lattice sets, at /3 = 7.192 and Nt = 24, corresponding to T = 1.5Tc. Also shown is the best 
fit to the form Eq. (|12|l with A = ST, and contributions of the different terms in the fit. LOG corresponds to the correlator 
constructed from the bcj"^ term in Eq. (|12|l and DIFF is the diffusive part of the correlator, constructed from the first term in 
Eq. (|12|l . (Right) The same information shown differently; the measured correlator, the best fit curve, and the diffusive part 
of the correlator are shown normalized to the leading order contribution. 



the plaquette value. This is translated to a temperature scale at = 8 using the information of (5c[Nr = 8) [2^ 
and two- loop scaling formula with a fitted correction (26j . Temperatures for other Nr arc easily calculated from the 
Nt ~ 8 temperature scale. The complete list of the lattice sizes, /3, and the corresponding temperature are shown in 
Table |TT] in the appendix, which also shows the parameters used in the multilevel algorithm for each (3. 



/3 6.76 6.80 6.90 7.192 7.255 
iV^ 20 20 20 24 20 
T/Tc 1.04 1.09 1.24 1.5 1.96 



TABLE I: List of lattices on which diffusion coefficients were extracted, and the temperatures. 



IV. RESULTS 



In order to calculate k, we calculated the electric field correlators, Eg. pil)) . for all the sets in table |TI1 From the 
correlators Ge{t), k can be calculated using Eq. ([9]). Use of the multilevel algorithm allowed us to get correlators at 
a few per cent level accuracy. In fact, we got < 2 — 3 % accuracy in all correlators except the two most central points 
of the 13 = 7.192, l.bT^ set. Fig. [5] shows Ge{t) for this data set. 

In order to get the momentum diffusion coefficient, k, we use the ansatz Eq. (|12[) for p{uj), and fit the Euclidean 
correlator using Eq. (|8]). It was not feasible to do a three parameter fit: the parameters A and k are strongly 
correlated. For a large range of A we can get very similar fit qualities. Instead, we fix A and get an estimate of k by 
doing a two-parameter fit. We discuss this further below and in Appendix IbI 

For the fit, minimization was carried out with the full covariancc matrix included in the definition of x^. We 
typically obtained acceptable fits to the correlators for ra in the range [Nt/A, Nt/2], with x^/d.o.f ^ 1. At shorter 
distances, lattice artifacts start contributing and the simple form of Eq. (fT2|) does not work well. Also using the 
leading order lattice correlator instead of the continuum form did not improve the quality of the fit. We, therefore, 
restrict ourselves to the long distance part of the correlator. 

In order to get a feel for the relative contributions of the different parts of the spectral function to the correlator, 
in Fig. [5] we show the correlators constructed from different parts of p{uj) separately. We take the best fit form of Eq. 
([T^ to the Nr = 24, 1.5Tc data set, for A = ST. The contributions to the total correlation function from the ui^ part 
of p{uj) and that from the diffusive part, the first term in Eq. (|12p . are calculated separately using Eq. (jH]). In Fig. 
[5] we have called these parts LOG and DIFF, respectively, and the correlator reconstructed from the fitted spectral 
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FIG. 3: (Left) The relative contribution of the diffusive part(DIFF) to the total correlator, compared to that of the leading 
order part (LOG), shown as a function of tT, for our different data sets. (Right) ois, defined through the scheme that the 
coefficient of the uj'^ term in p(u]) is 8as/9. 



function has been called Fit. The correlator is seen to be dominated by the contribution from the huj^ term over the 
virhole range of distance. Hovifever, the diffusion term has a substantial contribution near the center of the lattice. In 
Fig. [2] it contributes nearly 20 % at tT = 0.5. This is seen more clearly in the right hand panel of Fig. [2l where the 
total correlator, the best fit, and the diffusive part arc shown normalized by the leading order correlator. Since the 
relative contribution of the diffusive part falls rapidly at shorter distances, it is difficult to get reliable estimates of k, 
with the usual assorted tests like stability with small change in fit range, for our smaller lattices with Nr = 12 and 
16. So in what follows, wc quote fit results for k only for our finer lattices, with Nr = 20 and 24 (Table U). 

For these lattices, we obtained stable fits for the central part of the correlator, with x^/d.o.f ^ 1 in all cases. We 
did a fully correlated fit by including the inverse of the full covariance matrix in the function to be minimized, 
whenever such a function was well-behaved. That turned out to be the case in all sets except the one at the 
highest temperature, the /3 = 7.255 set in Table HI In this case we used an uncorrelated fit for our best estimate. The 
difference between the correlated and the uncorrelated fit was included in the systematic error. Our results for k/T'^ 
at various temperatures, using the ansatz Eq. (jl2p and A = 3T, arc shown in Fig. 21 The statistical error, shown by 
the solid (red) band, is obtained from a jackknife analysis. 

The choice of A = 3T for the central value was based on the fact that in all the sets, with A ~ 3T, the diffusion term 
contribution to the spectral function, aw in Eq. (jl2p . is numerically close to the large w term, bui^ , when the diffusion 
term sets in (i.e., at cj = A). Admittedly, this choice is somewhat arbitrary. In fact, the main source of uncertainty in 
our fit estimate, shown by the dashed (green) band in Fig. |4l is A. To estimate the possible error introduced through 
our central value of A = 3T, we varied A in the range [2T, oo). We also looked at the fit form Eq. ([T5)) . and did the 
same exercise with it. The details of the fit results for Eqs. (|12I13|) and various A are given in Appendix |B] Quite 
often the fit value for these variations comes outside the statistical error band of Fig. 21 A systematic uncertainty 
band is therefore introduced, of sufficient size so as to include the central fit values for all these variations. 

For the correlation functions of gluonic observables, major finite volume effects have been observed if the spatial 
size of the lattice is so small that some of the spatial directions get deconfined; on the other hand, at least for spatial 
correlation functions, finite size effects are small when the transverse directions are not deconfined [l^. To avoid large 
finite size effects, we choose lattices such that the spatial directions are confined. Since the electric field correlator 
also has contribution from the low uj part, it could be more sensitive to finite volume effects. However, as we discuss 
in appendix |B1 the correlation functions do not show any significant finite volume effect even when LT ^ 2. Therefore 
we do not expect large finite size corrections to our results obtained from lattices with LT > 2. 

It is instructive to look at the relative contribution of the diffusive part to the total correlator at different distances. 
In the left panel of Fig. [31 we show the correlator coming from the diffusive part of the fitted spectral function, 
normalized by the leading order part, for all the lattices of Table HI The notation is similar to that used in Fig. ^ 
except here we show the 1 a band and not the best fit value. At all temperatures, except the one at the highest 
temperature, the diffusive part is seen to reach about 5 % level by tT ~ 0.3. Note that the accuracy of our correlator 
is better than this. Also no significant trend of temperature dependence is seen in this figure. This is, of course, 
translated to the lack of significant temperature dependence of k/T^ in this temperature regime. Fig. [H 
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FIG. 4: The momentum diffusion coefficient, k, in units of T , shown as a function of temperature in the temperature range 
Tc < T < 2Tc. The error bars witli (red, solid) Une sliow the jackknifed error. The (green, dashed) error bars are an estimate 
of the size of the various systematic uncertainties, as discussed in the text. 



b, the coefficient of the term in p{u!) 
function is 



is also of some interest. In perturbation theory, the leading order spectral 



(15) 



To get an idea of the strength of the coupling at these temperatures, we use Eq. and the fit coefficient b, Eq. 
([T^ . for a nonperturbative estimate of a^. The estimate of as obtained this way is shown in Fig. [31 If Ge{t) is the 
properly normalized current, then the NLO calculation of Ref. [l^ can be used to connect this as to a^'^(/x). It is 
interesting to note that the coupling is rather small, about 1/4 near and going down to ~ 0.18 at 2 T^. This is 
in rough agreement with a similar measurement in Ref. [2l| from fit to vector current correlators, and other, more 
detailed, calculations of as at such temperatures from static observables [13] . 

In order to present our calculation in the context of RHIC, it seems convenient to use the Einstein relation between 
the diffusion coefficient, and k. 



D 



T 



Mt]d 



2T2 

K 



(16) 



In Eq. (|16|) rjD is the drag constant. In Fig. [5] we show the diffusion coefficient in the temperature range Tc '^T '2Tc, 
obtained using Eq. (|16p . The solid (red) error bar is the statistical error from a jackknife analysis. The bigger 
errorbars show the range of values covered by the different systematics analyzed in Table [B] 

Two points are worth noting in this figure. First, the nonperturbative value of the diffusion coefficient is rather 
small in the temperature range considered. In the next section we discuss in more detail the comparison with 
perturbation theory, but the diffusion coefficient shown here is nearly an order of magnitude smaller than the leading 
order perturbation theory. Second, there is no strong temperature dependence, at least in the temperature range 
Tc<T < 1.5Tc. 



DISCUSSION 



In this work, we studied the momentum diffusion coefficient, k, of heavy quarks in a gluonic plasma. As mentioned 
in the introduction, the large elliptic flow of the heavy flavor mesons, seen in the PHENIX experiment at RHIC, 
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FIG. 5; The diffusion coeiBcient, plotted as DT, in the temperature range ^ {Tc,2Tc). The error bars with solid red hues 
show the statistical error. The dashed(green) error bars are an estimate of the size of the various systematic uncertainties, as 
discussed in Appendix [B] 



seems to be well explained in a Langevin framework, if k, is large. Perturbation theory seems to be unstable for this 
quantity in the temperature regime of interest for RHIC physics, and the leading order PT prediction is at least an 
order of magnitude too small to explain the experimental data. Our aim in this work was to calculate the momentum 
diffusion coefficient nonperturbatively, and to see if the deviation from perturbation theory is of the size required to 
explain the experimental data. 

Using the formalism of Refs. [H, calculated k from the correlator of E'", the electric field operator. From 

the Matsubara correlator of the electric field, k was calculated through Eq. (j9]) using the ansatz for p{u!), Eq. (|12p . 
In order to compare our results with the perturbative calculation of 0] and experiments @, we used Eq. ([T5)) to get 
the diffusion constant, D. The diffusion coefficient so obtained is found to be considerably smaller than the LO PT 
estimate Q. For high temperatures such that itld/T <^ 1, the leading order estimate of DT is Q 



367r r A 2r 1 C'(2)\ Nf / 4T 1 C'(2)\ 

Cpg^ I \ mo 2 C(2) / 2 V f^D 2 C(2) / 

where Cp — {N^ — l)/2Nc is the color Casimir and Nf is the number of flavors of thermal quarks. At very high 
temperatures, DT diverges as ^/a^- As one comes down in temperature, Eq. (|17p is not reliable any more and one 
needs to use the complete leading order estimate. To get this, we use Eq. (11) of Q, with as(3T) determined using 
the plaquctte measurement js^l and mo taken from lattice measurements [27j . For example, at 1.5 Tc, ag^^{3T) w 
0.23 and mo/T w 2.345, leading to DT ~ 14. A similar calculation at 2.25 T^ and 3 Tc yield DT ~ 18.5 and 21, 
respectively, for the gluon plasma. A comparison with Fig. [S] reveals that this is almost an order of magnitude larger 
than the nonperturbative result for the gluon plasma. 

Interestingly, while Eq. ([TT)) seems to have a strong dependence on Nf, on putting values for the different quantities 
the Nf = 2 results are numerically not very different at similar values of T/Tc- The next-to-leading order (NLO) 
contribution to the diffusion constant has also been calculated in perturbation theory. At similar temperatures, with 
as ^ 0.2, this gives the DT - 8.4/(27r) for Nf = 3 ^. While a similar reduction for Nf = will bring the NLO 
PT result much closer to the nonperturbative estimate, it is rather disconcerting to find that the NLO result differs 
by almost an order of magnitude from the LO result. Indeed, one clearly will have to resort to calculations of higher 
orders/resummations before taking the perturbative estimates seriously. 

As already mentioned, there have been other attempts to calculate the diffusion coefficient using lattice gauge 
theory, so far only in the gluon plasma. In Ref. [l^ , preliminary results for an extraction of the diffusion coefficient 
from the vector current correlator c7jC was presented. The value of DT found at 1.5 T^ was considerably smaller 



(17) 
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FIG. 6; The difFusion coefficient of Fig. [5l shown here as 2-kDT, in the temperature range Tc < T < 2Tc. Also shown is the 
raiiee preferred by the V2 measured by PHENIX [^. The band is obtained from a comparison of Fig 40 of Ref. 0| and Fig 4 
of Q). The LO PT value at 1.5 Tc 0] is also shown. 



than LO PT, and is smaller than our results at that temperature, though consistent within systematics. Ref. [Ta| 
has also attempted extracting k/T^ from the electric field correlator, Eq. ([7]), but with considerably different analysis 
strategy. This calculation, which was concentrated mostly on considerably higher temperatures, found a very small 
value of k/T^, which docs not agree with ours in the temperatures where we overlap. On the other hand, a very 
recent calculation [l7| . which also focusses mostly at higher temperatures, is in much better agreement with ours in 
the temperature range of overlap. 

The heavy quark diffusion coefficient has also been calculated in a very different theory, the A/" = 4 supersymmetric 
Yang-Mills theory at large 'tHooft coupling XtH = ctgiVc, using AdS/CFT correspondence [H, [l^. In fact, a large 
part of the formalism used by us was introduced in Ref. For the A/" = 4 SYM theory for Nc oo and large XtH, 
Ref. [13 gets 

OT.'^(^)\ (18, 



2n \XtH 

Note that the dependence of D on the coupling in Eq. (|18p is parametrically different from that in Eq. (fT7|) . Of 
course, this theory is very different from QCD in many respect. Moreover, it exploits crucially symmetries which 
QCD docs not have. However, to get a feel for what kind of value such a functional dependence would give, one 
can somewhat arbitrarily put parameters relevant for QCD in Eq. (jl8p . Setting Nc = 3 and as = 0.23, one obtains 
DT ~ 0.2 from Eq. ([T8| . which is lower than, but in the same ballpark as our estimate. 

Our results are for quenched QCD, i.e., there are only thermal gluons but no thermal quarks in our fireball. So 
a comparison with experimental results needs to be done with care. A conservative approach would be to say that 
comparison of the results in Fig. [5] with the perturbative results for quenched QCD give us an indication of how 
much the nonperturbative results can change from the perturbative results in the deconfined plasma at moderate 
temperatures < 2Tc. Even then, the results are most encouraging since they indicate that the nonperturbative 
estimate for DT can easily be an order of magnitude lower than LO PT, bringing it tantalisingly close to values 
required to explain the V2 data. 

In a bit more optimistic fashion, one can hope that our results, as plotted in Fig. [Sj will be even quantitatively close 
to a similar figure in full QCD when it is computed. The reason for such a hope is that dimensionless ratios of various 
quantities arc known to scale nicely between quenched and full QCD if plotted as function of T/Tc- Also the LO PT 
result, Eq. (|17p . shows such a trend. In this spirit, in Fig. |6]we compare the lattice results with the experimental data. 
The lattice results seem to be a little above the best fit value for PHENIX, though reasonably close within our large 
systematics. Interestingly, our lattice results seem to show very little temperature dependence in the temperature 
regime studied here. For comparison, we also show in the same plot the leading order PT result for quenched QCD 
(see estimate below Eq. (|17p. which is, of course, very far from both the nonperturbative result and the experimental 
value. 

The most straightforward direction for possible refinement of our calculation is, of course, to go to finer and 
bigger lattices. A nonperturbative calculation of the rcnormalization constant will also be of great help in accurate 
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quantitative prediction. The nontrivial next step would be the inclusion of the light thermal quarks in the calculation. 
The multilevel algorithm cannot be directly used in that case, because of the nonlocality of the quark determinant. 
It would be an interesting challenge to come up with better ways to obtain similarly precise results even in the full 
QCD case. 



VI. ACKNOWLEDGEMENT 

We would like to thank Sourcndu Gupta for numerous discussions, and for insightful comments on the manuscript. 
The computations were done in the framework of Indian Lattice Gauge Theory Initiative (ILGTI) . The brood cluster 
of the department of theoretical physics, TIFR, and the gauge and chiral clusters of the department of theoretical 
physics, lACS, were used for this work. We would like to thank Ajay Salve and Kapil Ghadiali for technical support. 
PM would like to acknowledge Department of Science and Technology (DST) grant no. SR/S2/HEP/0035/2008 for 
the cluster chiral. The research of SD is partially supported by a Ramanujan fellowship from DST. The research of 
RVG is partially supported by a J. C. Bose fellowship from DST. 



[1] R. Baler, et al, Nucl. Phys.B 483 (1997) 291. 
[2] Y. Dokshitzer & D. Kharzeev, Phys. Lett.B 519 (2001) 199. 
[3] G. D. Moore and D. Teaney, Phys. Rev.C 71 (2005) 064904. 
[4] M. G. Mustafa, Phys. Rev.C 72 (2005) 014905. 

[5] A. Adare et al. (PHENIX CoUab.), Phys. Rev. Lett.98 (2007) 172301. 

B.I. Abelev et al. (STAR CoUab.), Phys. Rev. Lett.98 (2007) 192301. 
[6] A. Adare et al. (PHENIX CoUab. ) , arXiv: 100_5_.1627 
[7] B. Svetitsky, Phys. Rev.D 37 (198 8) 2484. 

[8] R. Rapp and H. van Hees. .arXiv:0903. 10961 in: R. C. Hwa (Ed.), Quark-Gluon Plasma 4, World Scientific. 

[9] S. Caron-Huot and G. Moore, J. H. E. P. 0802 (2008) 081. 
[10] See, e.g., H.B. Meyer, Eur.Phys.J. A47 (2011) 86, for a recent review. 
[11] P. Petreczky and D. Teaney, Phys . Rev.T) 73 (2006) 014508. 
[12] H-T. Ding, et al.. larXiv: 1107.0311 1 

[13] J. Casalderrey-Solana and D. Teaney, Phys. Rev.D 74 (2006) 085012. 

[14] S. Caron-Huot, M. Laine and G. D. Moore, J. H. E. P. 0904 (2009) 053. 

[15] Y. Burnier, M. Laine, J. Langelage and L. Mether, J. H. E. P. 1008 (2010) 094. 

[16] H. B. Meyer, New^.JPhys. 13 (2011) 035008. 

[17] A. Francis, et al.. [aFXi7:1 109.3941 

[18] J. Kapusta and C. Gale, Finite Temperature Field Theory, Cambridge University Press. 
[19] G. P. Lepage and P. Mackenzie, Phys. Rev.D 48 (1993) 2250. 
[20] Y. Koma and M. Koma, Nucl. Phys.B 769 (2007) 79. 
[21] H-T. Ding, et al., Phys. Rev.D 83 (2011) 034504. 

[22] M. Laine, G. Moore, O. Philipsen and M. Tassler, J. H. E. P. 0905 (2009) 014. 
[23] M. Liischer and P. Weisz, J. H. E. P. 0109 (2001) 010, J. H. E. P. 0207 (2002) 049. 
[24] S. Gupta, Phys. Rev.D 64 (2001) 034507. 
[25] G. Boyd, et al., Nucl. Phys.B 469 (1996) 419. 

[26] R. G. Edwards, U. M. Heller and T. R. Klassen, Nucl. Phys.B 518 (1998) 377. 

[27] O. Kaczmarek and F. Zantow, Phys. Rev.D 71 (2005) 114510. 

[28] S. S. Gubser, Nucl. Phys.B 790 (2008) 175. 

[29] S. Datta and S. Gupta, Phys. Lett.B 471 (2000) 382. 

[30] Note that the observed flow depends not only on the difi'usion coefiicient but also on various other details like the geometry 

of the collision, evolution, equation of state of the plasma, etc. See @] for a comprehensive review. 
[31] The efflciency is, of course, dependent on both /3 and Nt. 

[32] The non-pert urbative value of as at the inverse lattice spacing scale, = (3.4/a)exp(— 5/6) is obtained from the plaquette 
values [ly] , and then the 2-loop beta function is used to flow to the scale fi — 3T. 



Appendix A: List of lattices, and details of the algorithm 

Here we list the lattices used in our calculations, and parameters for the multilevel calculation. The parametrs 
for the multilevel calculation are also given. The last three columns correspond to the number of sublattices the 
lattice was divided in, the number of sublattice averaging between measurements (# update), and the number of 
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measurements of each correlation function (# conf). 
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TABLE II: Details of the lattices used in the calculation. Also given are the parameters for the multilevel algorithm for each 
set. 



Appendix B: Details of various systematics discussed in Sec. IIVI 

In this section we discuss in some detail some of the systematic uncertainties mentioned in Sec. IIVI 

• Uncertainties in the ansatz for the spectral function: 

In Sec. |ll]we have discussed the ansatz for spectral function used by us to get the diffusion coefficient. As we 
have discussed there, we have no first principle handle on the form, and have used forms for the low-w part 
motivated by other studies. In particular, we have introduced a cutoff in Eqs. ((T2l[T3)) . As we discussed in Sec. 
IIVI the quality of the fit is rather insensitive to A: as we vary A, we get a different best fit value for k, but for 
a range of A. the does not change appreciably. This is probably an example of the zero mode solutions we 
discussed in Sec. |TT1 We illustrate this in Fig. [71 for the {(3 = 7.192, 1.5Tc) set. In the figure the contribution of 
the diffusive part of the correlator is shown for the best fit parameters at various A. The corresponding values 
of K, shown in Table iBl vary substantially as A is varied; however, the total contribution of the diffusive part to 
the correlator hardly changes as we change A from 2T to 4T. 

Of course, the cutoffs in Eqs. ()12I13|) are an approximation: one does not expect discreet jumps in p(a;). It is 
not unreasonable to expect, however, that changing the sharp cutoff with a smooth one will not change things 
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FIG. 7: The change in the relative contribution of the diffusive part of the (3 = 7.192, 1.5 Tc data set, as we change A in Eq. 
(|12|l . The notation is same as Fig. [2] 



significantly and that the flat direction we encounter is of more general origin. For the purpose of this work, we 
take the conservative approach of letting A vary between [2r, oo), and include the values of k thus obtained in 
the systematic uncertainty band. We consider this range to be conservative because for A < 2T, plugging back 
the fit solution to construct p{uj), we get a large jump at cj = A, since at this value acu in Eq. (jl2p is much 
bigger than 6^"^. In order to quote a central value for the fits, we investigated for what value of A aw ~ buj^ 
for w = A. For all the sets of Table U this happens around A ^ 3T. Therefore, we use this value of A to quote 
the central value. Admittedly, this criterion is arbitrary, and the green band in Figs. 21 [5] is probably the more 
robust object. 

In Table IB] we also repeat this exercise of varying A for the fit form P2{^)- When the fit values obtained are 
outside the systematic band, the band is extended to include them. 
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TABLE III: Fit form dependence of k/T^. 



• Fit range and fit quality dependence: For the fit values quoted in Table |B1 we have used the range T,nin to 
Nr/2, where Tmin is the smallest r for which we got a good x^, and the is defined using the full covariance 
matrix. In all sets except one, we could get a good with Tmin > A'r/4, and increasing Tmin slightly did not 
change the fit value significantly. The set where we could not get such a stability with Tmin was the set at /3 = 
7.255. In this case only the uncorrelated i-e., the diagonal covariance matrix in the definition of x^, allowed 
such stability. So for this set, we used the uncorrelated fit value for our central estimate. The difference between 
the uncorrelated and the correlated best fits is then taken as an additional source of systematic uncertainty in 
this case. In fact, the systematic uncertainty band for this set in Fig. |4]is dominated by this contribution. 

• Finite volume effect: We explored finite volume effects by looking at LT ~ 2-4 on some of our coarser lattices, 
and smaller variations of LT in two of our finer lattices. At the correlation function level itself, no statistically 
significant finite volume effect could be seen once LT > 2. To make this statement quantitative, we do a x^ 
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comparison of the long distance part of the correlator, which should be the most sensitive to finite volume effects. 
For the correlator calculated on two lattices at the same /3 and but different Ng, we construct the quantity 



^-/4 V'^i(r)2+a2(r) 
For the different sets in Table this quantity is listed below. 
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Here the third column shows the LT values of the lattices whose correlators are being compared. At the level 
of accuracy of our correlators, we do not see any significant finite volume effect for LT = 2. So we believe our 
results, calculated on admittedly small lattices, will not be severely affected by finite volume effects. Even for 
the LT=1.8 set at /3 = 6.9, Nr = 20, where the correlator does show a statistically significant effect, Table [B] 
reveals that the error in k, due to finite volume effect is smaller than our other systematics. 

• Renormalization factor: 

The lattice correlator is multiplied by a renormalization factor, Eq. ([TT|) . to get Ge{t). Clearly, an error in 
Ze will affect k multiplicatively. As we discussed in Sec. in the absence of a nonperturbative evaluation of 
the renormalization factor, we have used the tadpole factor, which takes into account the quadratic self energy 
correction of the gluon lines, to renormalize G^^*. 

A perturbative renormalization factor, using heavy quark effective theory, has been calculated in Ref. (l7j . 
Below we tabulate the two renormalization factors for the different lattice spacings we have: 





6.76 6.80 6.9 7.192 7.255 


T/T, = 


1.04 1.09 1.24 1.5 1.96 


^tad 


1.230 1.232 1.226 1.210 1.207 


^HQET 


0.831 0.832 0.834 0.841 0.842 



Over the temperature range of interest to us, there is a near-constant factor 1.43-1.48 between the two renor- 
malization schemes. We do not include such a factor in our band of systematics, since it is easy to convert our 
results to Z^'^^'^ factor. We note that while this indicates a rather large reduction in k of order 30-32 %, it will 
not change our qualitative conclusions. 



